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1.  Introduction 


Since  their  introduction  in  the  late  1980’s,  essentially  non-oscillatory  (ENO)  and  weighted 
essentially  non-oscillatory  (WENO)  schemes  [1,  2,  3,  4,  5,  6]  have  become  two  of  the  standard 
methods  for  approximating  the  solution  to  hyperbolic  equations.  They  have  been  shown  to  be 
quite  successful  in  a  variety  of  applications.  ENO  and  WENO  schemes  approximate  the  solu¬ 
tion  by  a  constant  on  each  grid  element,  and  advance  in  time  explicitly  by  using  a  higher  order 
polynomial  reconstruction  of  the  current  solution.  The  ENO  scheme  computes  several  high  order 
reconstructions,  each  on  a  different  stencil,  but  uses  only  the  stencil  on  which  the  solution  is  the 
smoothest  (presumably  the  stencil  does  not  cross  a  shock).  The  WENO  schemes  exploit  the  fact 
that  often  a  convex  linear  combination  of  all  the  reconstructions  can  achieve  even  higher  order 
accuracy  than  ENO  at  a  specific  point  when  the  solution  is  smooth.  A  nonlinear  weighting  of  the 
polynomial  reconstructions,  based  on  the  solution  smoothness,  is  actually  used  to  avoid  stencils 
that  cross  shocks  in  the  solution. 

Because  reconstruction  stencils  need  to  avoid  shocks,  ENO  and  WENO  are  essentially  upwind 
methods,  and  some  form  of  flux- splitting  is  necessary  when  the  system  of  equations  is  subject  to 
flow  in  both  directions.  The  central  methods  avoid  flux  splitting,  and  central  WENO  (CWENO) 
schemes  were  developed  later  in  the  1990’s  [4,  6].  These  methods  use  staggered  grids,  offset  by 
one-half  grid  element,  and  so  require  higher  order  WENO  reconstructions  at  the  midpoint  of  a  grid 
element. 

Formally  fifth  and  even  ninth  order  CWENO  (CWEN05  and  CWEN09)  schemes  have  ap¬ 
peared  that  use  the  classic  WENO  methodology  to  define  high  order  reconstructions  from  low 
order  ones  [4,  6].  However,  there  is  no  classic  CWEN03  scheme  in  the  literature,  due  to  the  need 
for  high  order  reconstruction  at  the  center  of  the  grid  element.  The  linear  weights  required  do  not 
exist  in  the  third  (and  seventh)  order  case.  Indeed,  quoting  from  Qiu  and  Shu  [6,  p.  194],  “We 
remark  that  it  is  easy  to  verify  that,  when  r  is  an  even  number  (r  =  2,4,.. .),  there  are  no  linear 
weights  that  satisfy  [the  appropriate]  condition  ...  .  Hence  there  are  no  third-,  seventh-,  eleventh-, 
etc.,  order  central  WENO  reconstructions  for  the  point  values  with  these  choices  of  stencils.”  A 
classic  CWEN03  scheme  would  combine  two  linear  approximations,  each  requiring  the  solution 
on  two  grid  elements,  for  an  overall  stencil  of  three  grid  elements.  However,  the  two  linear  poly¬ 
nomials  agree  at  the  center  point,  and  so  no  convex  combination  can  be  higher  than  second  order 
in  a  classic  WENO  reconstruction. 

Third  order  CWENO  schemes  have  been  proposed  in  the  literature  by  Levy,  Puppo,  and 
Russo  [7,  8].  We  will  call  these  schemes  “augmented,”  because  they  use  a  convex  combination  of 
the  two  linear  approximations  augmented  directly  by  a  higher  order  polynomial  (a  quadratic)  in 
the  case  of  one  space  dimension.  Thus  there  is  no  weighting  of  low  order  polynomials  to  achieve 
higher  order  accuracy,  as  in  classic  WENO.  Rather,  high  order  accuracy  comes  directly  from  the 
quadratic  reconstruction.  One  might  view  these  methods  as  some  kind  of  adaptive  ENO  schemes, 
because  they  choose  the  linear  weights  somewhat  arbitrarily  and  use  the  WENO  smoothness  indi¬ 
cator  to  bias  the  stencil  combination  from  the  quadratic  in  smooth  regions  to  one  of  the  linears  in 
non-smooth  regions.  To  put  it  another  way,  classic  WENO  methodology  uses  lower  order  stencils, 
and  weights  them  to  increase  accuracy  in  regions  where  the  solution  is  smooth,  as  opposed  to  the 
augmented  CWEN03  schemes,  which  include  the  higher-order  stencil  with  the  low  order  stencils 
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and  weights  them  to  reduce  accuracy  near  discontinuities.  Clearly  both  approaches  have  a  great 
deal  of  merit. 

In  this  paper,  we  present  a  classic  CWEN03  scheme  that  indeed  combines  only  two  linear 
polynomials.  It  also  uses  a  compact  stencil  of  three  grid  element  values  in  the  reconstruction.  The 
key  to  this  work  is  to  develop  a  technique  that  allows  us  to  reconstruct  average  solution  values  on 
a  grid  different  from  the  original.  In  the  case  of  CWEN03,  two  natural  reconstruction  grids  are 
defined,  one  given  by  shifting  the  computational  grid  by  one-half  grid  element  and  resizing  by  a 
factor  of  one-quarter  locally,  and  the  other  is  nonuniform  and  described  later.  Re-averaging  (or 
re-mapping)  the  solution  from  the  original  computational  grid  to  the  new  reconstruction  grid  is 
then  accomplished  using  high  order  integration,  as  described  in  a  recent  paper  [9]  (see  also  [10, 
11]).  Once  this  re-averaging  to  the  new  grid  is  complete,  we  need  a  third  order  reconstruction 
of  the  solution  at  the  endpoints  of  a  reconstruction  grid  element,  rather  than  at  the  center  of  a 
computational  grid  element.  The  linear  weights  always  exist  for  the  endpoints  of  a  grid  element, 
and  so  we  can  proceed  to  define  a  classic  CWENO  scheme  that  is  formally  third  order  accurate. 
A  formally  seventh  order,  and  in  fact  any  odd  order,  scheme  could  be  devised  with  the  same 
technique,  although  we  do  not  present  the  details  here  (see,  e.g.,  [12,  13]  for  very  high  order 
schemes). 

When  linear  weights  exist  in  a  classic  WENO  reconstruction,  they  may  still  be  problematic 
if  some  of  them  are  negative  [11].  Fortunately,  Shi,  Hu,  and  Shu  provided  a  technique  to  handle 
this  case  [14].  Our  new  re-averaging  technique  provides  an  alternate  way  to  deal  with  negative 
weights,  by  avoiding  them  altogether.  One  can  re-average  to  a  shifted  and  perhaps  locally  rescaled 
reconstruction  grid  that  includes  the  targeted  point  as  an  endpoint,  since  the  linear  weights  for  the 
endpoint  case  are  known  to  be  positive. 

Standard  WENO  reconstructions  cannot  be  used  for  arbitrary  points,  since  we  may  not  have 
the  required  linear  weights.  However,  in  some  cases  one  may  wish  to  obtain  formally  high  order 
approximate  values  at  arbitrary  points,  say  in  a  postprocessing  step.  Our  re-averaging  technique 
can  be  used  to  overcome  this  limitation.  This  is  done  in  [15],  wherein  an  Eulerian-Lagrangian  (or 
semi-Lagrangian)  scheme  is  devised  that  requires  high-order  reconstructions  at  arbitrary  points. 

We  also  present  in  this  paper  an  extension  of  the  scheme  to  problems  in  more  than  one  space 
dimension.  It  turns  out  to  be  relatively  straightforward  to  do  so,  in  the  sense  that  much  of  the 
computation  remains  the  same  as  in  the  one  dimensional  case;  that  is,  the  scheme  is  developed  in 
an  essentially  tensor  product  form,  iterating  the  one  space  dimensional  scheme.  We  therefore  only 
present  the  case  of  two  space  dimensions.  Since  it  is  built  on  one-dimensional  computations,  the 
scheme  is  relatively  straightforward  to  implement,  and  it  is  computationally  efficient,  as  opposed 
to  the  multi-dimensional  approach  of  Levy,  Puppo,  and  Russo  [7,  8].  We  also  prove  that  the 
scheme  is  indeed  formally  third  order  accurate. 

The  re-averaging  WENO  reconstruction  technique  is  presented  in  the  next  section.  Section  3 
presents  the  new  CWEN03  scheme  for  a  scalar  conservation  equation  in  one  space  dimension. 
The  scheme  is  extended  to  problems  posed  in  two  space  dimensions  in  Section  4.  Numerical 
results  in  one  space  dimension  are  presented  in  Section  5,  and  the  technique  is  applied  to  the  Euler 
system  in  Section  6.  Numerical  results  in  two  space  dimensions  are  presented  in  Section  7,  and 
the  paper  ends  with  conclusions  in  the  final  section. 
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Figure  1:  The  re-averaging  process.  The  three  values  i<;,  and  Mi+i  are  re-averaged  to  v,  which  is  a  third  order 
approximation  to  the  average  of  u  on  the  subinterval  [x,_  1/2. ^]- 


2.  The  re-averaging  technique 

Let  us  fix  the  notation  for  the  computational  grid  by  choosing  the  spacing  h  >  0  and  setting  the 
grid  points  to  be  x,-_i/2  =  (i  -  1  )h  and  the  central  points  to  be  x,  =  (i  -  1/2 )h.  The  z'th  grid  element 
is  then  ((/  -  1  )h,  ih \  =  [*,-1/2,  *1+1/2].  and  its  center  is  x,.  The  function  u(x)  is  approximated  on  the 
grid  by  its  element  averages 

|  nxi+ 1/2 

Ui  ~  -  I  u(x)dx.  (1) 

h  Jx- 1/2 

Our  re-averaging  WENO  reconstruction  technique  is  a  two  stage  process.  In  the  first  stage, 
we  define  the  reconstruction  grid  and  re-average  the  computational  grid  element  averages  m,  to  our 
new,  reconstruction  grid.  We  will  take  a  reconstruction  grid  that  is  a  shift  and  local  rescaling  of  the 
computational  grid.  The  re-averaging  is  done  according  to  the  high-order  integration  as  described 
in  [9].  The  second  stage  of  the  process  is  to  apply  any  standard  WENO  reconstruction  for  which 
the  linear  weights  exist  for  the  new  reconstruction  grid. 

In  this  section  we  discuss  the  first  stage  of  the  process.  The  second  stage  is  determined  by 
the  application,  and  is  described  in  the  next  section  when  we  define  our  new  CWEN03  schemes. 
More  details  of  the  re-averaging  can  be  found  in  [9,  10,  11,  15],  but  we  present  here  only  the  two 
natural  cases  needed  for  a  CWEN03  scheme.  To  fix  ideas,  suppose  that  we  need  to  reconstruct, 
say,  the  function  value  at  x„  the  center  of  grid  element  i,  to  third  order  accuracy  using  only  the 
values  M;_i,  Uj ,  and  ui+  \ .  Again,  there  is  no  standard  third-order  WENO  reconstruction  in  this  case, 
since  the  linear  weights  do  not  exist  for  the  center  point. 


2.1.  The  re-averaging  process 

The  key  to  the  overall  WENO  reconstruction  process  is  to  produce  a  high  order  reconstruction 
of  the  function  averages  on  a  new  grid.  Such  a  high  order  reconstruction  was  described  for  fifth 
order  reconstructions  in  a  recent  paper  [9].  We  present  the  details  for  third  order  reconstruction 
here. 

As  illustrated  in  Fig.  1,  consider  a  point  £  such  that  x,-_  1 /2  <  <$  <  *,+i/2,  which  gives  a  subinterval 
[Xi- 1/2,  <$]  of  the  central  element  [ .v  11  /2,  *r+i/2]-  Our  goal  is  to  reconstruct  u  using  only  iz,_i ,  ii,.  and 
ui+ 1  so  that  we  obtain  a  high  order  approximation  to  the  integral 

Ui(x):=  I  u(y)dy. 

dxi- 1/2 


The  ENO  approximations  are  linear  functions  [3,  11],  one  on  a  left  stencil  using  1  and  un 
giving 

„  _  Xi  -  X  _  X  -  X  _1 

PL(x)  :=  w,_  1 — - —  +  Uj - . 


h 
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and  one  on  a  right  stencil  using  w,  and  «,-+1,  giving 


_  xM  -x  X-Xi 

Pr(x)  :=  Uf - - - +  ui+ 1- 


h 


h 


We  need  to  integrate  these  to  determine  the  linear  weights.  However,  it  is  simpler  to  determine 
the  integrals  by  approximating  directly  the  primitive  function  C//(a).  The  integrals  of  PL  and 
PR  are  simply  the  two  quadratic  polynomials  interpolating  U,  using  the  proper  masses.  The  left 
polynomial  approximates  at  a/_3/2,  0  at  a,_i/2,  and  /?«,■  at  a/+i/2,  and  it  is  given  by 


„  ,  .  _  (X~  Xi-1/2)(X  -  Xi+1/2)  _  (X-  Xi- 3/2)(a  -  Xi- 1/2) 

Ql(x)  :  =  —  M|_  j - — - +  Ui - 


2/? 


2/7 


-f 

*JXj- 


PL(y)dy. 


1/2 


Similarly,  the  right  polynomial  approximates  0  at  jcx-_i/2,  hit,  at  a,-+i/2,  and  //(w,  +  m;+1)  at  a,-+3/2.  It  is 


_  (a  -  A/_i/2)(a  -  A/+3/2)  _  Ax-Xi-  1/2)(a  -  A/+i/2) 

2/?W  :=  -77/ - 7 -  +  (77/  +  77/+i) - 


h 


2/7 


-f 

JXi- 


PR(y)  dy. 


1/2 


These  two  quadratics  are  combined  using  the  linear  weights  yL  and  yR  =  1  -  yL  into  the 
polynomial 

Q(x)  :=  yLQAx )  +  7r27?(a). 

To  set  yy,  we  compare  2(a)  to  the  higher  order  cubic  polynomial  C( a)  that  interpolates  all  the 
data,  i.e.,  it  is  -hu^\  at  a ,-_3/2,  0  at  a,-_i/2,  /7m;  at  a,-+i/2,  and  h(Uj  +  77,+1)  at  a;+3/2.  Thus 


C( a)  :=  77 


(a  -  A/_i/2)(a  -  A/+i/2)(a  -  A/+3/2)  _  (a  -  a,_3/2)(a  -  a,_i/2)(a  -  A/+3/2) 


7-1  ■ 


+  (77/  +  77/+ 1) 


6h2 

(a  -  A/_3/2)(a  -  A/_i/2)(a  -  A/+i/2) 
6/72 


2/72 


If  we  wanted  to  approximate  the  primitive  function  to  high  order  at  a  =  we  would  choose  yy 
so  that  g(£)  =  C(£),  no  matter  the  values  of  77,_i,  m,,  and  m/+1.  This  overdetermined  system  has  a 
solution,  and  the  resulting  linear  weights  are  positive,  and  given  by 

7  At)  =  >  0,  7fM)=  1  -  Jl(A)  =  >  0,  p(^):=^~f~1/2  e[0,l].  (2) 

3  3/7 

These  linear  weights  ensure  that  the  formal  order  of  accuracy  of  <2(£)  to  C//(£)  is  0(hA).  That  is, 
using  the  ratio  of  distance  across  the  zth  grid  element  p  =  piA),  we  have 
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2  “,+1J 

h. 


(3) 


In  WENO  applications,  the  linear  weights  yL(£)  =  (2  -  p)/3  and  y«(^)  =  (1  +  p)/3  would  be 
modified  by  the  smoothness  indicator,  as  described  below. 

However,  we  want  to  approximate  an  average  of  u  between,  say,  =  a,-_i/2  +  p\h  and  - 
A/- 1/2  +  PiK  which  is 


v(pi,p2) 
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Figure  2:  The  uniform  reconstruction  grid.  The  computational  grid  is  labeled  by  x,  while  the  reconstruction  grid  is 
typeset  below  it  and  labeled  in  The  three  values  u,  \ ,  S,-,  and  ui+\  are  re-averaged  to  v_i,  vq,  Vi,  and  v2>  which  in 
turn  form  two  sets  of  three  values  for  third  order  WENO  reconstruction  at  x,  =  ft. 


This  involves  differences  of  the  primitive  function,  and  the  linear  weights  do  not  always  exist  in  this 
case  [9].  However,  for  our  CWEN03  scheme,  we  need  averages  of  u  only  over  certain  subintervals 
of  the  zth  grid  element  for  which  the  linear  weights  exist,  and,  in  fact,  they  are  positive.  We  present 
the  two  natural  cases  of  reconstruction  grids  in  the  next  two  subsections. 

2.2.  Case  1,  a  uniform  reconstruction  grid 

We  take  in  the  first  case  of  a  new  CWEN03  scheme  the  reconstruction  grid  shifted  by  one- 
half  grid  element  and  resized  by  a  factor  of  one-quarter  locally.  That  is,  as  depicted  in  Fig.  2,  the 
new  grid  near  the  target  point  xt  covers  the  interval  [a',_i/2,  xi+ 1/2],  the  zth  element,  and  is  given  by 
ft 2  =  Xi-1/2,  fti  =  xi- 1/2  +  h/ 4,  ft  =  xi  =  *z- 1/2  +  h/ 2,  ft  =  v,-+1/2  +  3/z/4,  and  ft  =  x,-+i/2. 

We  need  to  re-average  the  solution  on  the  computational  grid  to  the  reconstruction  grid  {ft2, 
fti ,  ft,  ft ,  ft}.  For  example,  to  compute  a  higher  order  approximation  to  the  function  average  over 
the  second  quarter  of  the  zth  grid  element,  i.e.,  to  get  v0  =  v(l/4, 1/2),  we  would  set  yL  and  yR  so 
that 


h 

-v0  sb  yi\Qifxi- 1/2  +  h/2)  -  Qifxi- 1/2  +  h/ 4)]  +  yz?[2z?(x/_i/2  +  h/ 2)  -  0/?(x,_i/2  +  /z/4)] 
=  C(Xj-]/2  +  h/2)  -  C(v,_i/2  +  /-z/4). 


The  result  is  that  for  the  quarter  subintervals, 
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Figure  3:  The  nonuniform  reconstruction  grid.  The  computational  grid  is  labeled  by  x,  while  the  reconstruction  grid  is 
typeset  below  it  and  labeled  in  The  three  values  w,_i,  u,.  and  S,+i  are  re-averaged  to  v_i  =  u,_\ ,  v0,  vi,  and  v2  -  Uj+i, 
which  in  turn  form  two  sets  of  three  values  for  third  order  WENO  reconstruction  at  x,-  =  £q. 


2.3.  Case  2,  a  nonuniform  reconstruction  grid 

In  this  second  case,  we  select  a  reconstruction  grid  that  is  nonuniform.  As  depicted  in  Fig.  3, 
we  cover  the  interval  [x,_3/2,  x/+3/2]  with  the  new  grid  f_2  =  x(_3/2,  =  x,_l/2,  =  xh  fx  =  xi+i/2 

and  =  x,-+3/2;  that  is,  we  simply  divide  the  center  element  in  half.  In  this  case,  formula  (3) 
applies  directly,  and  leads  to 
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2.4.  Modification  of  the  linear  weights  for  non-smoothness 

Of  course  to  avoid  discontinuities  within  the  WENO  framework,  we  modify  these  linear 
weights  with  the  usual  smoothness  indicator 


ISl 


if . . 

C=  1  Uxi- 1/2 


r  2(-l 


d'Pfx) 

dx{ 


i2 


dx,  k  =  L,R, 


(13) 


using  the  original  linear  polynomials.  For  a  uniform  grid,  an  explicit  integration  yields 

IS'L(u )  =  ( Uj  -  Ui- 1)2  and  ISR(u)  =  (ui+ x  -  Hi)2. 

The  nonlinear  weights  are  then  given  by 


(14) 


7l 


fflL, u)  := 


(e  +  /Si(n))2 


7 l 


+ 


1  -  7  L 


and  y‘R(yR,  u)  :=  l yR ,  n).  (15) 


(6  +  /5:(n))2  (6  +  /5'(n))2 


where  e  >  0  is  some  small  parameter.  That  is,  in  (5)— (12),  the  linear  weights  in  the  formulas  are 
replaced  by  the  nonlinear  weights  y/(yL ,  u)  and  Yr(jr,  it),  using  the  original  linear  weights  yL  and 
yR  as  given  on  the  right-hand  side. 

In  (5)— (1 1),  we  have  the  re-averaged  values  of  v,  defined  in  (4),  for  the  subintervals  that  we 
will  need  later,  accurate  to  0(h3)  when  the  function  u(x)  is  smooth,  since  b_  -  (i  is  0(h).  Note  that 
the  re-averaging  uses  only  the  linear  polynomials  Pl(x)  and  PR(x)  (in  integrated  form  QL(x)  and 
Qr(x)),  and  that  the  stencil  is  compact,  using  only  the  grid  element  averages  m(_,  ,  uh  and  uM.  It  is 
also  essentially  non-oscillatory  near  shocks,  since  it  drops  back  to  only  using  one  linear  stencil  by 
the  nonlinear  smoothness  indicator  correction  to  the  linear  weights. 


Remark.  For  the  quarter  subintervals,  the  re-averaged  function  values  (5)-(8)  do  not  preserve 
mass  when  the  nonlinear  weights  are  used,  since  the  linear  weights  differ  between  subintervals.  A 
mass  conservative  re-averaging  is  possible  [15]  if  one  computes  the  integrals  always  from  the  left 
endpoint  x,_i/2.  That  is,  as 
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wherein  the  linear  weights  have  been  converted  to  the  nonlinear  weights  y'L,  y'R  using  (14)— (15). 
Since  the  nonlinear  weight  corrections  are  done  consistently  for  each  subinterval,  mass  is  con¬ 
served.  However,  mass  conservation  at  this  stage  is  not  required  for  our  CWEN03  scheme,  since 
the  re-averaged  function  values  are  used  merely  to  obtain  high-order  approximations  to  the  func¬ 
tion  at  a  point.  Of  course,  the  overall  scheme  itself  does  conserve  mass.  In  computations,  we  did 
not  see  any  noticeable  difference  between  using  (5)— (8)  and  (16)— (19).  Of  course,  the  half  subin¬ 
terval  formulas  use  the  same  weights,  and  so  conserve  mass  at  this  and  all  stages  of  the  scheme. 

Remark.  Both  of  our  two  re-averaging  procedures  exhibit  a  superconvergence  effect  for  the  higher 
order  reconstruction  of  a  point  value.  To  be  precise,  if  one  re-averages  [m,_i,  uh  «/+,  }  to  the  uni¬ 
form  or  nonuniform  reconstruction  grid  to  obtain  [v_i,  Vo,  Vi,  v2}  and  then  reconstructs  u(xj)  in  the 
standard  way  either  from  the  left  stencil  {v_i,  v0>  vi}  or  the  right  stencil  {vo,  Vi,  v2]  using  the  linear 
weights  directly  without  nonlinear  modification  in  both  stages  of  the  reconstruction,  then  one  ob¬ 
tains  the  accuracy  0(h4).  To  prove  this  fact,  one  needs  to  show  that  the  procedure  is  exact  when  the 
underlying  function  u(x)  is  a  cubic  polynomial.  In  fact,  it  is  enough  to  show  that  when  u(x)  =  x3. 


8 


Figure  4:  An  illustration  of  a  time  step  in  a  CWENO  scheme.  The  shifted  grid  element  average  2  at  the  advanced 
time  t,H  1  is  computed  from  the  reconstructed  solution  K,(x)  and  /?,+i (x)  at  time  t",  which  gives  the  initial  mass  of  the 
solution  in  place  over  the  region  [x,,  x,+ j],  and  the  flux  values  over  time  at  the  grid  locations  x,  and  xi+  \ . 

the  reconstruction  value  is  zero.  Conceptually,  when  the  linear  weights  are  used,  the  reconstruc¬ 
tion  value  is  the  same  as  what  arises  from  the  quadratic  approximation  of  u(x )  evaluated  at  x  =  0 
in  both  stages  of  the  reconstruction.  For  the  re-averaging  stage,  since  the  cubic  is  an  odd  function 
and  the  reconstruction  grid  is  symmetric,  the  quadratic  approximation  of  x 3  must  also  be  odd,  i.e., 
a  multiple  of  x.  The  second  stage  will  reconstruct  the  value  zero  for  x,  and  the  proof  is  complete. 
Of  course,  one  can  simply  verify  the  result  for  u(x)  =  x3,  by  applying  the  formulas  (5)-(8)  and 
(26)-(27)  for  the  uniform  reconstruction  grid  and  (9)— (12)  and  (28)-(29)  for  the  nonuniform  re¬ 
construction  grid.  This  superconvergence  result  will  break  down  when  nonlinear  weights  are  used. 
Nevertheless,  one  can  check  the  accuracy  of  a  computer  code  implementation  of  the  algorithm  by 
checking  for  this  superconvergence. 


3.  CWEN03  for  a  scalar  conservation  law  in  one  space  dimension 

Given  /'(«;  x,  t ),  consider  the  initial  value  problem  for  a  hyperbolic  conservation  law 


du  df(u ) 

—  +  =  0, 

dt  dx 

u(x,  0)  =  uq(x), 


t  >  0, 

xgI. 


(20) 

(21) 


A  description  of  the  details  of  a  Central  WENO  scheme  can  be  found  in  [4,  6].  We  present 
herein  a  summary  of  the  entire  scheme  for  completeness,  although  the  only  new  feature  is  the  way 
we  handle  the  WENO  reconstruction  at  the  center  of  the  grid  element  by  using  our  re-averaging 
procedure. 

Let  t°  =  0  <  tl  <  t2  <  ■  ■  ■  define  the  time  levels.  In  a  CWENO  scheme,  the  grid  shifts  one-half 
grid  element  each  time  step.  We  describe  a  step  from  time  tn  to  tn+ 1  using  the  grid  indexing  we 
have  been  using,  i.e.,  xt  is  the  center  of  an  element  of  the  grid  at  time  t".  CWENO  is  based  on 
integrating  (20)  over  a  shifted  grid  element  [xh  xi+i]  in  the  space  variable  x  and  [f,  tn+1]  in  time  t. 
That  is,  as  illustrated  in  Fig.  4,  the  solution  average  on  the  shifted  element  [jc,-,  X/+ \  ]  at  time  tn+l  is 
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m"++',  and  it  is  set  from  the  approximation 


-n+ 1 


u 


(+1/2 


w(a,  tn)  dx 


£  fML,dt+ 


(22) 


The  first  integral  in  (22)  uses  the  WENO  reconstruction  function  RJx,  t")  of  {m,_i  ,  uh  w/+1,  w,+2}. 
To  be  more  precise,  Ru(x,  f )  is  defined  by  R\(x)  on  [x/^1/2,  x(-+i/2],  which  involves  1,  if,-,  w,+i  },  and 
7?,+i(.v)  on  [.V/+1/2,  a',+3/2J,  which  involves  {«,,  w;+i,  m/+2}.  The  required  integrals  were  actually  given 
above  in  (1 1)  and  (10),  and  they  are 
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(23) 


wherein  the  linear  weights  happen  to  all  be  1  /2,  but  the  two  halves  of  the  integral  use  a  different 
position  for  the  smoothness  indicator  to  define  the  nonlinear  weights. 

For  the  two  flux  integrations,  i.e.,  the  second  and  third  terms  in  (22),  we  use  a  Gaussian  quadra¬ 
ture  rule  for  the  integration  of  the  flux  in  time  over  [/",  t'H  1  ].  and  the  two-point,  third  order  rule  is 
appropriate  here.  Let  A  f  =  t"+l  -  f  and 


1 

2 


2V3 


(24) 


For  j  =  i,i  +  1,  the  quadrature  rule  is 

r'”+1  a  tn 

J  m\xj  dt  *  —[mxj,  f  +  A/"0g))  +  fUdXj,  f  +  Afe+G))\.  (25) 

Therefore,  the  main  task  remaining  is  to  find  properly  reconstructed  values  for  u(xj,  tn  +  A/"0=), 
i.e.,  at  the  centers  of  the  computational  grid  elements  in  space  and  at  the  Gauss  points  in  time. 


3.1.  Flux  evaluation  in  space  using  re-averaging 

Evaluation  of  (25)  in  space  uses  our  new  re-averaging  technique.  To  obtain  high  order  accuracy 
for  u"  ~  u(xj,  /"),  we  can  use  either  the  uniform  or  nonuniform  grid  from  Sections  2.2  and  2.3.  We 
re-average  the  solution  {u"_v  m”,  uni+l)  to  the  chosen  reconstruction  grid  to  obtain  {v_i,v0,vi,v2}. 
As  illustrated  in  Figures  2  and  3,  there  are  actually  two  possible  stencils  of  three  re-averaged 
quantities  for  standard  WEN03  reconstruction  to  the  point  xr  That  is,  we  can  define  two  linear 
WENO  reconstruction  polynomials,  one  for  each  of  the  stencils  {v_! ,  v0,  iq }  and  {v0,  iq ,  v2}. 

For  the  uniform  reconstruction  grid,  we  have  the  left  and  right  stencil  reconstructions 


<L:=r«(  1/3,  v) 
<*:=ri(: 2/3,  V) 


-  jv-i  + 
1  1 
r° + 2"1 


+  r^(2/3,v) 
+  7r(  1/3,  v) 


L2Po  +  2Pl 


:Vi 


-v2 


(26) 

(27) 
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using  the  smoothness  indicator  defined  at  the  appropriate  index  for  differences  of  values  of  v.  For 


the  nonuniform  reconstruction  grid,  it  is  easy  to  derive  the  linear  weights,  and  we  have 

uIl  :=  Tl(1/4^)  -  ^v_i  +  ^v0  +  7r(3/4,  v)  ^v0  +  ^vi  ,  (28) 

uur  '■=  rl(3/4,  v)  ^v0  +  ^vi  +  Tr(1/4,  v)  ^vi  -  ^v2  •  (29) 

In  this  latter  case,  it  should  be  noted  that  the  smoothness  indicator  must  be  computed  from  (13) 
using  our  nonuniform  grid.  The  left  stencil  has 

IS'L{v)  =  |~(v0  -  v_i)j  and  IS‘R(v)  =  (vj  -  v0)2,  (30) 

and  for  the  right  stencil, 

ISlL( v)  =  (vj  -  v0)2  and  ISlR(v)  =  |^(v2  -  vi)j  .  (31) 

To  produce  a  central  scheme,  we  average  the  results  of  these  two  stencils  and  define  the  formally 
third  order  reconstructions 

u(xh  f)  *  <  :=  and  /(m(x„  r»);  x„  r")  *  .//'  :=  /(«”;  x,-,  r").  (32) 


(As  an  alternative,  one  could  introduce  flux-splitting  and  take  the  upstream  weighted  stencil,  but 
we  present  only  the  central  scheme  herein.) 


3.2.  Flux  evaluation  in  time  using  Runge-Kutta 

Evaluation  of  (25)  in  time  is  handled  in  the  usual  way,  in  our  case,  by  using  a  two-stage  Runge- 
Kutta  method  combined  with  the  natural  continuous  extension  of  Zennaro  [16]  applied  for  fixed 
x  =  x,  to  the  initial  value  problem 


du  df(u ) 


and  u(xi,F)  =  uni. 


One  needs  to  approximate  the  space  derivative  of  /'(«;  x,  t)  using  the  reconstructed  point  values  f". 
The  first  stage  of  the  Runge-Kutta  method  is  to  define 


g(i)  :=  -fL(ii2,r)fLJk  ~yR(U2,r) 


-.ir ,  am 

h  dx 


(so  that  the  overall  scheme  has  a  compact  stencil  of  five  points)  and  then  set  the  intermediate  values 

wt  :=  ul  +  Atng{l).  (34) 

The  second  stage  is  then  to  define  f"+]  =  f(w,\  x,-,  t"+l)  and 
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’  h 


dfiu) 

dx  Xi,tn+l 


The  Runge-Kutta  method  would  be  completed  by  defining  m”+1  :=  u "  +  — -(g(1)  +  g(2)).  However, 

we  do  not  need  u  at  t  =  tn+l  ,  but  rather  at  the  Gauss  points  in  time.  The  natural  continuous 
extension  of  the  Runge-Kutta  values  over  the  time  interval  gives  us 


u{xuf  +  A  fOf)  *  u';+e±G  :=  a'l  +  ^[(2  Of  -  (0§)2)g(1)  +  (0§)2  g(2)], 
5.5.  Summary  of  the  scheme 

Combining  the  steps  above,  (22)  with  (23)  and  (25)  simplifies  to 


(36) 


utli2-.=  fL(\ll,un) 


- 1%  i  +  tM! 
8  1  8  ' 


+  yf\i/2,un) 


1  3 

~U-  +  — 

8  1  8  1 


+  7r(1/2,  un) 

+  %+1(l/2,  un) 


- u "  +  -u" 

8  '  8  1 


-  1. <"  +  A<"e0) + /(«";, "G,* 


2h 

At"-\ f(un+<l 

2  h 


+  —  [/(« 


G-  x„  r"  +  ArX)  +  /(«r"G;  T,  f  +  Afe+C)\, 


(37) 


where  G  are  defined  in  (33)-(36),  which  requires  u"  and  fj  defined  in  (26)-(32)  and  v  defined 
either  in  (5)-(8)  for  a  uniform  reconstruction  grid  and  in  (9)— (12)  for  a  nonuniform  reconstruction 
grid.  One  also  needs  Of-  defined  in  (24)  and  and  y'K  defined  in  (14)— (15). 


Remark.  The  nonuniform  reconstruction  grid  is  somewhat  more  computationally  efficient.  How¬ 
ever,  the  numerical  results  (see  Example  5.2)  perhaps  indicate  that  the  uniform  reconstruction  grid 
is  the  better  choice. 


4.  CWEN03  for  a  scalar  conservation  law  in  higher  space  dimensions 


We  consider  now  extension  to  problems  with  more  than  one  space  dimension.  Extension  to 
three  and  even  higher  dimensions  is  straightforward  from  the  case  of  two  space  dimensions,  so  we 
present  only  this  case,  which  is 


du  df(u)  dg(u ) 

dt  dx  dy 

u(x,y,0)  =  u0(x,y). 


(x,y)  €  M2,  t  >  0, 
(x,y)  6  M2, 


(38) 

(39) 


for  the  given  f(u\ x,y,  t )  and  g(u;  x,y,  t). 

The  computational  grid  is  a  tensor  product  of  uniform  grids.  The  spacing  in  x  is  denoted  hx  >  0 
and  the  spacing  in  y  is  hy  >  0,  and  we  tacitly  assume  that  the  ratio  hx/hy  is  bounded  above  and 
below  as  these  parameters  are  refined.  Again,  for  a  step  from  f  to  tn+1,  the  center  points  of  the 
grid  elements  are  (x,-,  yf  =  (. ihx,jhy )  and  the  current  approximation  to  the  element  averages  are 
denoted 

|  nyj+ 1/2  nXi+i/2 

~  rr  u(x,y,  f)  dxdy. 

hXhy  Jyj.1,2  Jxi-i/2 
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(40) 


Integration  of  the  differential  equation  leads  us  to  the  basis  of  our  approximation  scheme,  to  wit 


-73+1 


1+1/2J+1/2  ~  h 


1  i  ryj+i  rXi+i 

—  I  u(x,y,tn)dxdy 

<chy  l  Jy  Jx 


x"y  v  Uyj  Uxi 

nyj+i 

m\ 

n-^+i 

g(u)  I 

. 


rf"+l  ryj+i  { 

x=XMdydt  +  J  J  f(u)\x=xdydt 


’ yj 

V+1  rxi+i 


dx  dt  + 


nr  nxj+i 

Z(Ut=y/ 

Jtn  Jx;  1 


dxdt). 


(41) 


4.1.  High  order  integration  in  two-dimensions 

The  first  double  integral  in  (41)  is  approximated  to  high  order  using  a  two-dimensional  ana¬ 
logue  of  (3).  To  derive  this  analogue,  let  x,_ i /2  <  d  <  xMp_  and  v,-_  1/2  <  //  <  v/+i/o.  Apply  (3)  to 
a  given  function  u(x,  y)  in  the  variable  x  keeping  y  fixed,  and  then  apply  it  again  in  y.  In  terms  of 
the  scaled  values 


£  -  Xj- 1/2 

h\ 


e  [0, 1] 


and  cr{rf) 


d  ~  yj-i/2 
hy 


6  [0,  1], 


(42) 


the  result  is 


■.=  r  f 

y 7 — 1/2  d 

jlM,  n) 


u(x,  y)  dx  dy 

yj- 1/2  oxi- 1/2 

P(  1  -P)cr(l  -  cr) 


p(l  -  p)  o-(l  +  cr)  _ 

+  - X - X - Ui-\J 


+ 


p(  1  +p)  cr(l  —  cr) 


M; /-I  + 


p(l  +  p)  cr(l  +  cr) 


+  7m(£  >1) 


+  7 rAZ,  d) 


+  7 r.r(^  d) 


P(  1  -P)  cr(3  -cr)_ 


1,7 


+ 


2  2  ",-w 
p(l  +p)  cr(3  -  cr) 


2  2 
P( 3  -p)cr(l  -cr) 


- 


p(l  -p)  cr(l  —  cr)  _ 

2  2  Ui~1J+1 
P(  1  +P)  cr(l  -  cr)  _ 

o  o  ui.j+] 


2  2 
p(l  -p)  cr(l  cr) 


P(3  -  p)  cr(  1  +  cr)  _ 
uUj- 1  +  o  o  uUj 


p(l -p)cr(l +cr)_ 
m<+i,/-i  o  o  M!+1J 


2  2  ^  2  2 
p(3  -  p)  cr(3  -  cr)  _  p(3  -p)  cr(l  -  cr)_ 

^ /j+l 


2  2 
p(l  -p)  cr(3  -  cr) 


Hj  ~ 


ui+l,j  + 


2  2 
p(l  -p)  cr(l  -  cr). 


2  2 

wherein  the  linear  weights  are  given  by  (2)  and 

7*,/(£  P)  :=  7*(£)  7c(A  k,£  =  L,  R. 


/zA,  (43) 


(44) 
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We  claim  that  the  accuracy  of  the  approximation  in  (43)  is  0(h5),  where  h  =  ma x(hx,hy). 
Let  us  write  (3)  in  terms  of  the  linear  approximation  operator  Az  and  the  error  Ez  -  0(h4)  when 
considering  an  integral  in  z  =  x  or  y,  e.g.,  for  jc  as 


u(x)  dx  =  Ax(u )  +  Ex(u)  . 


Let  Uj(y)  =  j-  fA,+l/2  u(x,y)dx  and  uAx) 

"t  JXi- 1/2  J 


kSyiZU{X'y)dy-  ThCn 


u(x,  y )  dx  dy 


[Ax(u(y))  +  Ex(u(y))\  dy 


-  Ax(  I"  u(y)  dyj  +  0(h5) 

'  dy.i-t/2  ' 

=  Ax(Ay(u )  +  Ey(u))  +  O(IE) 

=  AxAy{u)  +  Ax{Ey{u))  +  0(h5), 


and 


Ax{Ey(u))  =  f  L'YiX.v))  J.r  -  Ex(Ey(u(x))  =  0(/z5), 

dxi-l/2 


since  the  error  term  is  smooth  when  u  itself  is  smooth,  and  so  the  claim  is  established. 

The  linear  weights  should  be  modified  by  the  two-dimensional  smoothness  indicator,  which 
can  be  written  in  terms  of  the  multi-index  a  as 


j  nyj+ 1/2  rw+i/2 

10-1=1,2  Jyj-i/2  Jxj 


IS!i  := 


ltd1' 


1  h2ai' 
ny 


1  [DiyPk/x,  y)f  dx  dy,  kJ  =  L,  R, 


(45) 


-1/2 


using  the  tensor  product  linear,  i.e.,  bilinear  polynomials  Pkj(x,  y)  that  match  the  data.  For  exam¬ 
ple, 


Pydx,  y)  :  =  H,-_ i. 


(Xi  -  x)(y  j  -  y) 


j-y 


+  ui,  j- 1 


h  xhy 


+  u 


i-i  J" 


(xi  -  x)(y  -  yh i) 


h  Ji, 


(x  -  Xi-i)(yj  -  y)  _  (x  -  Xi-i)(y  -  yh 0 


h xhy 


+  M 


hj 


hxhy 


An  explicit  integration  yields 

7 

75 ^z(a)  =  (, Uij  -  Ui-ij)1  +  (iiij  -  1)2  +  —  (M/,J  -  M/-IJ  -  +  M/-1J-1)2,  (46) 

and  the  other  smoothness  indicators  75  ^(w),  75  ZfZ  («).  75 ^(w)  can  be  computed  by  symmetry.  The 
nonlinear  weights  are  then  given  by  rescaling  as  in  (15). 
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We  approximate  the  first  integral  in  (41)  according  to  the  decomposition  into  four  grid  elements 


ryj+i  nxi+i 
Jyj  Jxi 


hxhy  Mi+ !  /2J+ 1  /2 


u(x,y,  t")  dxdy 


pyj+i  nxi+i 

dyi+ind  Xj 


u(x,y,  f)  dxdy  + 


r%yj+ 1/2  nxi+ 1 
Jv;  */*/. 


)V+ 1/2  *SXi+i/2 
~*yj+ 1  nxi+x/2 


u(x,y,  f )  dx  dy 


+ 


pVj+l  M 
dyj+1/2  d  Xj 


>yj  ^Xi+i/2 

"'J’j+l/2  nXi+1/2 


u(x,y,  f)  dxdy  + 


oyj+1/2  /’jc 

Vi 


«(x,  y,  f")  <7.v  Jy  (47) 


by  applying  (43)  with  p  =  cr  =  1  /2,  or  its  symmetric  equivalent,  to  these  four  integrals  individually. 
In  terms  of  the  generic  approximate  integral  operator  (divided  by  hxhy) 

U (i,  j,  u ;  i'i,  12,13;  71,72,73)  := 

_  r  1  3_  3_  9_ 

7^(1/4,m)^— il/ij,  +  —  uiuh  +  —uhJl  +  —  uhj2 

~ij  r\  ia  -\\  ^  -  1-  15  _  3  _ 

+  74(1/4,  u)y—Uiuh  -  —uiuj3  +  —  ui2j2  -  —ui2j3 

A  5  15  1  3 

+  74/ !/4, «)[—«,,,/,  +  —uhj2  -  —uhJl  -  —uhj2 

~ijn,A-J 25-  5  _  5  _  1_ 

+  74(1/4,  u)  —uhj2  -  —uhJ 3  -  —  uh,h  +  —Ui 3J3 


(48) 


the  result  is 


M'j+\ /2j+]  12  :=  +1,7  +  1,  w" ;  U  i  +  1,  i  +  2 ;  7, 7  +  1,7  +  2) 

+  f/(i  +  1,7,  w"  ;  i  +  1,  i  +  2 ;  7  +  1,7, 7  -  1) 

+  C/(i,  7  +  1,  m"  ;  i  +  1,  /,  i  -  1 ;  7, 7  +  1,7, 7  +  2) 

+  17 (i,  7,  w'1 ;  i  +  1,  i,  i  -  1 ;  7  +  1, 7, 7  -  1).  (49) 

4.2.  Reconstruction  of  the  solution  at  the  current  time  tn  at  the  center  points 

In  anticipation  of  our  need  to  evaluate  the  flux  terms  in  (41)  in  space  and  time,  we  first  re¬ 
construct  to  high  order  the  solution  at  the  center  points  of  the  grid  elements  at  the  current  time. 
We  do  this  by  first  re-averaging  and  reconstructing  in  y  and  then  doing  the  same  in  x.  Of  course, 
we  could  do  this  the  other  way  around,  and  we  could  even  compute  using  both  approaches  and 
average  the  results  to  preserve  symmetry.  However,  for  computational  efficiency,  we  simply  use 
the  stated  approach. 

For  each  index  i,  re-averaging  in  y  produces  the  four  values  {v,-  _i ,  vi;0,  v,-,i ,  v,^},  using  either  (5)- 
(8)  for  the  uniform  reconstruction  grid  or  (9)— (12)  for  the  nonuniform  reconstruction  grid,  but  of 
course,  using  the  nonlinear  weights.  These  re-averaged  values  are  then  reconstructed  using  (26)- 
(27)  or  (28)-(29)  on  two  stencils  to  produce  the  two  values  v^L  and  which  are  averaged  to 
obtain  v(  J.  These  reconstructed  values  are  high  order  accurate  approximations  of  the  averages  in  v, 
i.e.,  of  y  «(v,y;,  /")  dx  (because  the  reconstruction  is  a  linear  operation  which  commutes 
with  integration  in  the  other  variable,  see  [9]  for  a  proof). 
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Now  that  we  have  high  order  accurate  averages  in  x  for  y  =  v7,  we  can  repeat  the  procedure 
in  .r,  that  is,  re-average  {v(_,  j,  v,-j,  vi+i j)  in  .v  using  (5)-(8)  or  (9)-(12),  reconstruct  using  (26)-(27) 
or  (28)-(29),  and  average  the  results  from  the  two  stencils  to  obtain  the  high-order  accurate  point 
value 

u'lj  ~  u(Xi,yhtn). 

Moreover,  we  define 


fu  ■=  Wp  xi>yj,  O  and  gl  -  :=  g(j/!.\  xhyj,  tn). 


We  remark  that  a  genuine  two-dimensional  reconstruction  of  the  center  values  could  be  devel¬ 
oped.  It  would  require,  say,  sixteen  two-dimensional  re-averagings  to  a  uniform  4x4  reconstruc¬ 
tion  grid  within  the  (/,  j)-grid  element,  and  then  a  two-dimensional  reconstruction  on  four  stencils. 
This  would  entail  much  more  work  than  using  the  approach  we  have  presented,  which  involves 
eight  one-dimensional  re-averagings  and  four  one-dimensional  reconstructions  (and  twice  this  for 
a  symmetric  implementation).  To  be  fair,  the  computation  of  (47)  would  be  trivial  in  the  former 
case,  and  requires  four  two-dimensional  re-averagings  in  the  latter  case,  but  the  work  is  still  less 
using  our  iterated  one-dimensional  approach. 


4.3.  Evaluation  of  the  solution  in  time  at  the  center  points 

To  handle  time  integration  of  the  flux  terms  in  (41),  we  again  use  a  two-point  Gauss  rule,  which 
requires  approximations  of  the  solution  «(jq,  yj,  t )  at  the  Gauss  times  t  =  f  +  A t"6^..  Following  a 
two-dimensional  analogue  of  the  usual  Runge-Kutta  procedure  with  the  natural  continuous  exten¬ 
sion  outlined  in  Section  3.2,  i.e.,  with  an  additional  term  involving  the  derivative  of  g  with  respect 
to  y,  we  obtain  our  approximations 

w  u(xhyj,f  +  A f%). 


and  set 


rn+8t. 

f  ■ 

J‘J 


fiun^-xuyhr% 


and 


n+8% 


8ij 


8(\ 


n+6 ; 


G;xi,yj,tn+ec). 


4.4.  Evaluation  of  the  flux  integrals  over  space  and  time 

The  evaluation  of  the  flux  integrals  in  (41)  is  accomplished  by  using  a  two-point  Gauss  rule  in 
time  and  a  WENO  reconstruction  of  the  space  integral.  For  example,  the  last  flux  integral  in  (41) 
is 


]  nt"+1  f'Xi+i 

vd  J,  s<u)U' 


dx  dt 


At" 


2/i  \h\ 


r*i+ i 

g(u)  I 

J  X; 


y=yj,t=tn+AtnG~ 


dx  + 


nxi+i 

g(u)\ 

Jxi 


y=yj,t=tn+Atn9^ 


dx 


We  would  like  to  use  (3)  for  the  approximation  of  the  space  integral,  but  we  are  using  high  order 
point  values  rather  than  element  averages  at  this  stage  of  the  computation.  Instead,  we  use  a 
standard  WENO  reconstruction.  The  space  integral  is  approximated  by  breaking  it  into  two  pieces 
at  .r  =  ;q+i/2.  For  example,  the  left  piece  over  [v,,  v,+  i/2]  is  reconstructed  using  {gi_ffgil  6 ,  g,  h , } 
by  matching  yL  times  the  integral  of  the  left  linear  interpolant  with  yR  times  the  integral  of  the 
right  linear  interpolant  to  the  integral  of  the  quadratic  interpolant.  The  linear  weights  turn  out  to 
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be  are  1/3  and  2/3.  Thus,  using  the  nonlinear  weights,  we  have  the  approximate  values  for  each 
time  level  as  <D yi+ip  i(g"+6£),  where 


Cia/s)  =  rTO1/3’  8(-)j) 


h 


+  g#U 


+  ri+1(2/3,g(.)j) 


1  3 

g<?<J  +  g^'+ij 


+  7«(2/3,  g(-)j) 
+  yR+\l/3,g,)J) 


1 

8' 


+  o  <?'+!.;' 


L  g  ■?<+ 1 J  g  <?<+2,/ 


Similarly,  define 


C+i/2(/)  =  ^-iri(!/3,  At))|  -  gib-i  +  g/ij  +7r(2/3 ,Ao) 


8‘ 


8“ 


+  yf1(2/3,A(,) 


1  3 

8  Ji’j  +  g/ij'+l 


+  7r+1(1/3,  Ac-)) 


3  1 

g/ij  +  g/id+1 
g  •/',./+ 1  —  gii'j+2 


(50) 


(51) 


Then  the  sum  of  the  four  flux  integrals  in  (41)  is 


<D 


n 

i+ 1.2, ;+ 1/2 


-yj  -  ^+u+i/2(r+^)  -  Ku+wir**)  + 

-  Kinj+ibr*)  -  Kuu+ibr1*) + ,2j(8n+e° 


(T^)  +  %+l/2(fn+d+o) 

)  +  %ipJgn+e+o)}. 


(52) 


4.5.  Summary  of  the  scheme  in  two  space  dimensions 

We  finally  summarize  the  scheme  based  on  (41).  We  first  use  (49)  and  (48)  to  compute  the 
mass  M"+] /2  j+\ /i  at  Ome  We  then  define  the  center  point  values  «/.,  /7,  and  g" j  using  the  re¬ 
averaging  and  reconstruction  procedure  of  Section  4.2,  and  then  produce  u"^.''0,  f”.  °G ,  and  g'.'  "0 
using  the  Runge-Kutta  with  natural  continuous  extension  procedure  of  Section  4.3.  Finally,  we 
compute  the  flux  0/+|  2  j+\/2  using  (50)— (52),  and  set 


-  n+ 1 


U 


1+1/2J+1/2  '■ 


Mn  +  cb" 

lvli+ll2,j+l/2  +  vy/+l/2J+l/2- 


(53) 


5.  Numerical  results  for  problems  in  one  space  dimension 

In  most  works,  the  nonlinear  weights  (15)  use  the  parameter  e  -  IE-6  [3].  However,  in  [8],  it 
was  found  that  augmented  CWEN03  works  better  when  using  larger  values  of  e  (say,  e  =  IE-4 
or  £  =  IE-2).  The  larger  the  value  of  e,  the  more  the  nonlinear  weights  resemble  the  linear  ones. 
Herein  we  sometimes  vary  e  from  its  usual  value.  In  the  second  example,  we  will  also  experiment 
with  changing  the  strength  of  the  measure  of  smoothness. 

All  our  test  results  use  periodic  boundary  conditions. 

5.1.  Example  1,  constant  linear  transport 

We  first  test  our  scheme  in  the  simple  case  of  constant  linear  transport  as  applied  to 

u,  +  ux  =  0,  x  e  (0, 2),  u0(x)  =  0.75  +  0.25  sin(7r.r). 
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Table  1:  Ex.  1,  linear  transport.  Error  and  convergence  order  at  T  =  10  with  At  =  0.2 /z  and  h  -  2/m.  We  use  e  =  IE-2 
in  the  nonlinear  weights  and  the  uniform  reconstruction  grid. 


111 

L\  error 

order 

error 

order 

10 

2.31287E-01 

1.78861E-01 

20 

4.67840E-02 

2.30560 

4.17714E-02 

2.09825 

40 

4.16134E-03 

3.49090 

4.72872E-03 

3.14299 

80 

4.89365E-04 

3.08806 

4.57875E-04 

3.36843 

160 

6.10702E-05 

3.00237 

5.04137E-05 

3.18307 

320 

7.63440E-06 

2.99988 

6.07320E-06 

3.05329 

640 

9.54460E-07 

2.99976 

7.52046E-07 

3.01356 

1280 

1.19299E-07 

3.00010 

9.37727E-08 

3.00358 

Table  2:  Ex.  1,  linear  transport.  Error  and  convergence  order  at  T  =  10  with  At  =  0.2 h  and  h  -  2/m.  We  use  e  =  IE-2 
in  the  nonlinear  weights  and  the  nonuniform  reconstruction  grid. 


m 

L'h  error 

order 

Lf  error 

order 

10 

2.26608E-01 

1.75287E-01 

20 

4.41577E-02 

2.35946 

3.92850E-02 

2.15767 

40 

4.11885E-03 

3.42235 

4.49674E-03 

3.12703 

80 

4.89164E-04 

3.07385 

4.47272E-04 

3.32966 

160 

6.10693E-05 

3.00180 

5.00641E-05 

3.15930 

320 

7.63440E-06 

2.99986 

6.06220E-06 

3.04586 

640 

9.54460E-07 

2.99976 

7.51702E-07 

3.01161 

1280 

1.19299E-07 

3.00010 

9.37620E-08 

3.00309 

Table  3:  Ex.  1,  linear  transport.  Error  and  convergence  order  at  T  =  10  with  At  =  0.2 h  and  h  =  2/m.  We  use  e  =  IE-4 
in  the  nonlinear  weights  and  the  uniform  reconstruction  grid. 


m 

Llh  error 

order 

L?°  error 

order 

10 

2.92394E-01 

2.25889E-01 

20 

1.64297E-01 

0.83161 

1.41641E-01 

0.67338 

40 

3.13130E-02 

2.39147 

3.56153E-02 

1.99167 

80 

3.99845E-03 

2.96925 

5.04667E-03 

2.81909 

160 

1.93416E-04 

4.36966 

2.84825E-04 

4.14718 

320 

9.51802E-06 

4.34490 

1.37065E-05 

4.37714 

640 

9.54468E-07 

3.31789 

9.91066E-07 

3.78974 

1280 

1.19299E-07 

3.00011 

1.01243E-07 

3.29116 

For  this  smooth  problem,  we  can  simply  use  the  linear  weights  (i.e.,  we  can  make  no  nonlinear 
smoothness  correction  to  the  linear  weights).  In  this  case,  we  observe  a  clean  third  order  rate 
of  convergence  for  the  scheme,  as  expected.  However,  the  scheme  will  normally  be  applied  to 
problems  that  have  non-smooth  behavior,  and  so  we  need  to  test  the  performance  of  the  nonlinear 
weights. 
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The  numerical  results  using  the  nonlinear  weights  are  not  as  clean.  We  show  L]h  and  Lf°  errors 
for  the  nonlinear  weights  with  a  large  e  =  IE-2  in  Tables  1-2,  comparing  the  uniform  and  nonuni¬ 
form  reconstruction  grids.  In  this  case,  we  see  a  fairly  clean  third  order  of  convergence.  We  also 
note  that  there  is  little  difference  between  the  results  for  uniform  and  nonuniform  reconstruction 
grids.  In  Table  3  we  show  the  same  comparison  but  using  a  smaller  e  =  IE-4  for  the  uniform  re¬ 
construction  grid  (results  for  the  nonuniform  grid  are  very  similar).  The  convergence  is  not  clean, 
but  it  appears  to  be  third  order.  Moreover,  the  errors  are  comparable  to  the  previous  case  of  a 
larger  e. 


Table  4:  Ex.  1,  linear  transport  with  the  alternate  initial  condition  uq(x)  =  sin4(7o).  Error  and  convergence  order  at 
T  =  1  with  At  =  0.2/;  and  h  —  2/m.  We  use  e  =  IE-2  in  the  nonlinear  weights  and  the  uniform  reconstruction  grid. 


111 

Llh  error 

order 

error 

order 

10 

6.51455E-01 

4.19944E-01 

20 

3.19967E-01 

1.02574 

3.64131E-01 

0.20574 

40 

8.84790E-02 

1.85452 

1.32501E-01 

1.45845 

80 

1.07544E-02 

3.04041 

2.11501E-02 

2.64727 

160 

9.36958E-04 

3.52080 

1.67371E-03 

3.65955 

320 

1.01509E-04 

3.20638 

1.27499E-04 

3.71449 

640 

1.24687E-05 

3.02522 

1.29981E-05 

3.29411 

1280 

1.55225E-06 

3.00588 

1.52999E-06 

3.08671 

We  also  test  our  scheme  with  the  initial  condition  Uq(x)  =  sin4f/r.v).  The  results  are  shown 
in  Table  4  (for  the  uniform  reconstruction  grid  and  e  =  IE-2).  The  scheme  gives  a  good  approxi¬ 
mation  in  this  case,  which  has  an  initial  discontinuity  in  the  derivative  at  x  =  1  and  x  -2. 

Finally,  we  remark  that  the  augmented  CWEN03  scheme  of  Levy,  Puppo,  and  Russo  [7,  8] 
show  results  similar  to  those  obtained  by  our  CWEN03  scheme  on  these  tests. 

5.2.  Example  2,  Shu ’s  linear  test 

We  take  a  standard  test  problem,  called  Shu’s  linear  test,  defined  for  x  6  [0,2],  and  solve  it 
for  various  choices  of  parameters.  The  main  results  use  a  mesh  of  m  =  1000  grid  points,  which  is 
reasonable  for  comparison  to  other  papers  using  CWEN05  with  about  200  points.  Our  CWEN03 
scheme  drops  to  second  order  near  the  discontinuities,  for  an  error  <9((2/1000)2)  =  0(4  x  10-6) 
versus  CWEN05,  which  is  third  order,  for  an  accuracy  of  0(( 2/200)3)  =  <9(10~6).  Results  are 
shown  after  linear  transport  by  speed  one  at  time  T  =  2.0  using  At  =  OAh. 

As  we  saw  in  the  previous  example,  a  larger  value  of  e  in  the  nonlinear  weights  (15)  can 
produce  a  cleaner  convergence  rate  on  a  smooth  problem,  since  this  biases  the  nonlinear  weights 
toward  being  the  linear  weights.  However,  a  larger  e  does  not  work  well  for  non  smooth  examples, 
such  as  Shu’s  linear  test,  since  the  use  of  the  unmodified  linear  weights  leads  to  oscillations  in  the 
solution.  As  an  alternative,  we  can  increase  the  strength  of  the  measure  of  smoothness. 

Instead  of  using  the  L2-measure  of  smoothness  as  given  above  in  (13),  we  use  the  //-measure 
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m  =  1000,  p  =  2,  nonuniform  m  =  1000,  p  -  4,  nonuniform 


Figure  5:  Ex.  2,  Shu’s  Linear  test.  The  solution  using  m  =  1000  grid  points  (h  =  0.002)  and  At  =  0.4 h  at  time  T  =  2.0. 
The  top  row  of  graphs  show  results  using  the  uniform  reconstruction  grid,  and  the  bottom  row  uses  the  nonuniform 
reconstruction  grid.  The  left  column  uses  p  -  2  in  the  nonlinear  weights,  and  the  right  column  uses  p  =  4.  The  results 
suggest  that  the  uniform  reconstruction  grid  and  the  greater  p  is  more  accurate. 


for  p  >  2.  That  is,  we  define 


if . . 

[=\  JXi-Up 


hpt~l 


d[Pk(x ) 
.  dxe 


p 


dx. 


k  =  L,R. 


(54) 


In  our  case  of  linear  polynomials,  the  new  measure  of  smoothness  is  given  simply  by  replacing 
the  exponent  2  by  p  in  (14)  for  the  uniform  grid  and  in  (30)  and  (31)  for  the  nonuniform  grid.  The 
nonlinear  weights  are  still  defined  by  (15),  using  expressions  like  (e  +  ISf')2,  k  =  L,R,  as  opposed 
to  the  more  standard  way  of  increasing  the  strength,  which  involves  the  expressions  (e  +  IS2'‘)P 
(see  [3]).  The  reason  for  using  this  alternate  way  of  increasing  the  strength  is  that,  in  our  case, 
increasing  p  does  not  change  the  strength  of  the  parameter  e. 

Our  main  results  compare  the  choice  of  p  in  the  nonlinear  weights  (using  a  fixed  e  =  IE-8)  and 
the  choice  of  reconstruction  grid.  As  shown  in  Fig.  5,  at  least  for  this  test  problem,  we  see  better 
accuracy  when  p  =  4  as  opposed  to  p  =  2.  We  also  see  better  accuracy  using  the  uniform,  versus 
the  nonuniform,  reconstruction  grid. 

We  also  varied  m  to  show  that  accurate  results  can  be  obtained  for  this  example  using  a  coarse 
grid,  as  long  as  p  in  the  nonlinear  weights  is  taken  to  be  large  enough  (p  =  6  or  p  =  8).  The 
results  appear  in  Fig.  6.  It  seems  to  be  an  open  question  as  to  the  best  choice  of  nonlinear  weight 
parameters.  The  case  p  =  2  is  generally  considered  to  be  adequate  [3],  but  it  is  not  the  best 
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m  =  800,  p  —  4 


m  =  800,  p  =  6 


m  =  400,  p  —  8 


Figure  6:  Ex  2,  Shu’s  Linear  test.  We  obtain  good  results  on  relatively  coarse  grids  provided  a  larger  p  is  used  in  the 
nonlinear  weights.  Here  the  uniform  reconstruction  grid  is  used. 


choice  in  this  example.  However,  the  Buckley-Leverett  test  of  Section  5.4  showed  no  significant 
difference  between  using  p  =  2  versus  p  -  4,  so  our  results  are  far  from  conclusive  in  the  general 
case. 

5.3.  Example  3,  Burgers’  equation 

In  this  example  we  test  Burgers’  equation  with  a  simple  initial  condition, 
ut  +  (w2/ 2)x  =  0,  x  €  (0, 2),  uq(x)  =  0.75  +  0.25  sin(^x). 

With  e  =  IE-4  in  the  nonlinear  weights,  we  see  third  order  convergence  in  Table  5. 


Table  5:  Ex.  3,  Burgers’  equation.  Error  and  convergence  order  at  T  =  1  with  A t  =  0.2 h  and  h  =  2/m.  We  use 
e  =  IE-4  in  the  nonlinear  weights  and  the  uniform  reconstruction  grid. 


m 

Llh  error 

order 

L“  error 

order 

10 

9.30494E-02 

1.22347E-01 

20 

2.48902E-02 

1.90242 

6.09936E-02 

1.00425 

40 

5.09583E-03 

2.28819 

2.82353E-02 

1.11116 

80 

9.462 10E-04 

2.42909 

8.57124E-03 

1.71992 

160 

1.56259E-04 

2.59822 

2.00149E-03 

2.09843 

320 

2.18843E-05 

2.83597 

3.53081E-04 

2.50300 

640 

2.82913E-06 

2.95146 

4.90373E-05 

2.84805 

1280 

3.52475E-07 

3.00477 

6.19008E-06 

2.98585 

We  now  study  Burgers’  equation  with  a  combination  of  a  shock  and  a  rarefaction.  The  initial 
condition  is  changed  to 


u0(x)  = 


0.5, 

1.0, 


if  x  <  0.3  and  x  >  0.75, 
if  0.3  <  x  <  0.75. 


With  only  m  =  80  grid  points,  we  see  good  accuracy  in  Fig.  7  compared  to  the  augmented 
CWEN03  [7,  8], 
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1.1- 


0.9- 
0.8] 

0.7- 

0.6- 

0.5- 

°-a.2  1  6  1  '0:2'  1  0:4'  1  W  1  1  o:s  1  1  1  '1:2 

Figure  7:  Ex.  3,  Burgers’  equation  with  a  shock  and  rarefaction.  The  solution  using  m  =  80  for  our  CWEN03 
scheme,  shown  in  black  squares,  and  the  augmented  CWEN03  scheme,  shown  in  blue  crosses.  (We  use  the  uniform 
reconstruction  grid.)  The  two  schemes  are  comparable  in  accuracy. 


’ot8~i  1  r“i 


t  —  0.1,  shock  and 
rarefaction  formation 


0.8- 

\ 

0.6 

'' 

0.6 

0.4- 

0.4- 

0.2 

0.2- 

t  =  0.2,  the  initial  pulse  t  =  0.4,  the  two  pulses  t  =  0.5,  the  two  pulses 
catches  the  advanced  pulse  merge  are  fully  merged 


Figure  8:  Ex.  4,  Buckley-Leverett  equation.  An  interaction  of  shocks  and  rarefactions,  resulting  from  the  evolution  of 
the  initial  condition  of  two  pulses  given  in  (55).  The  red  line  is  the  reference  solution,  given  by  CWEN05  with  a  very 
small  h  -  1/1280  and  At  =  1/15360.  The  black  squares  are  our  CWEN03  results  using  m  -  80  and  At  =  0.2 h. 


5.4.  Example  4,  Buckley-Leverett  equation 

This  last  example  involves  (20)  with  the  Buckley-Leverett  flux  function 


/(«)  = 


u2  +  (1  -  u)2 


and  an  interaction  of  shocks  and  rarefactions.  The  initial  condition  is 

( 1  -  20x  for  0  <  jc  <  0.05, 

“oW  =  1 0.5  for  0.25  <  x  <  0.4,  (55) 

lo  otherwise, 

and  has  two  pulses  that  merge  in  time.  We  use  m  =  80  grid  points  and  p  =  2  in  the  nonlinear 
weights  (no  major  difference  was  observed  using  p  =  4).  The  results  are  shown  in  Fig.  8,  and  the 
scheme  handles  the  merging  of  the  two  pulses  quite  well  for  using  only  80  grid  points. 
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6.  Application  to  the  Euler  system 


For  a  polytropic  gas,  the  energy  is  E  -  p/(y  -  1)  +  pu1  / 2,  where  p,  p,  and  u  are  the  particle 
pressure,  density,  and  velocity,  and  y  is  the  adiabatic  index  (y  =  (/+2)//  =  1.4,  where  /  =  5  is  the 
number  of  degrees  of  freedom  of  each  gas  particle).  The  one-dimensional  dynamics  is  described 
by  the  Euler  equations 


P 

pu 

E 


+ 


pu 

pu2  +  p 
u{E  +  p) 


(56) 


To  respect  the  characteristic  structure  of  the  system,  it  is  usual  to  reconstruct  the  needed  infor¬ 
mation  using  the  linearized  system.  We  expand  (56)  into  a  system  of  the  form  Ut  +  A(U)UX  =  0, 
that  is,  into 


p 

pu 

+ 

l  E  J 

t 

\ 

0 


1 


0 


\(y  ~  3 )u2  (3  -  y)u  y  -  1 

u  i  \  3  u(E  +  p)  ( E  +  p )  , ,  2 

,(y  -  1)m - - (y  -  1  )u-  yu 

P  p 

The  linearized  equation  can  be  diagonalized  as  A0  =  L°  A(U°)  (L0)-1,  giving 

1  m°  -  c°  0  0 


'  P  A 
pu 


+ 


0 

0 


u°  0 
0  u°  -  c° 


A  (  p  \ 

pu 
E 


=  0, 


P 

pu 

E 


=  0. 


where  c  =  sj(yp)/p  is  the  sound  speed.  The  re-averaging  of  Section  2  and  the  reconstructions 
needed  for  the  integral  over  the  previous  time  level,  i.e.,  (23),  are  performed  on  the  diagonalized 
variables  L°(p,pu,E).  We  then  transform  back  to  the  original  variables  (p,pu,E)  by  applying 
(L0)-1,  and  the  time  evolution  of  the  scheme  continues  for  the  nonlinear  system  (56),  as  described 
above  in  Section  3. 


6.L  Example  5,  Riemann  problems  for  the  Euler  equations 

For  this  series  of  tests,  we  specify  a  discontinuous  initial  condition,  written  in  terms  of  the 
primitive  variables  p,  u,  and  p.  As  is  typical,  we  only  report  the  density  p;  the  other  variables  show 
comparable  accuracy. 

Sod’s  1-D  shock  tube  test.  The  one-dimensional  shock  tube  test  of  Sod  uses  the  initial  condition 

[pi  =  l,ui  =  0,pi  =  1,  for  x  <  1/2, 

p,  u,p  =  < 

\pr  =  l/S,ur  =  0,pr=l/10,  for  x  >  1/2. 

The  results  are  shown  in  Fig.  9.  Good  results  are  obtained  using  m  =  200  and  m  =  400  grid  points. 


Lax’s  1-D  shock  tube  test.  The  one-dimensional  shock  tube  test  of  Fax  uses  the  initial  condition 

{pi  =  0.445,  ut  =  0.698,  pi  =  3.528,  for  x  <  1  /2, 
p,  u,p  =  < 

|pr  =  0.5,  ur  =  0,  p,  =  0.571,  for  x  >  1/2. 

Reasonably  good  results  are  shown  in  Fig.  10,  using  m  =  200  and  m  =  400  grid  points. 
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Figure  9:  Ex.  5,  Sod’s  1-D  shock  tube  test.  The  density  profile  at  time  T  =  0.16  using  At  =  0.1  h  and  m  =  200  (left) 
and  m  =  400  (right)  grid  points. 


Figure  10:  Ex.  5,  Lax’s  1-D  shock  tube  test.  The  density  profile  at  time  T  =  0.16  using  At  =  0.05/;  and  m  =  200  (left) 
and  m  =  400  (right)  grid  points. 


Figure  11:  Ex.  5,  Arora  and  Roe’s  Mach  3  shock  tube  test.  The  density  profile  at  time  T  =  0.09  using  At  =  0.05/;  and 
m  =  200  (left)  and  m  =  400  (right)  grid  points. 


Arora  and  Roe ’s  Mach  3  shock  tube  test.  The  mach  3  shock  tube  test  of  Arora  and  Roe  has  the 
initial  condition 


[pi  =  3.857,  uj  =  0.92, pi  =  10.333,  for  x  <  1/2, 
p,u,p  =  < 

[pr  =  1,  ur  =  3.55,  p,  =  1,  forx>  1/2. 

Good  results  are  shown  in  Fig.  11,  using  m  =  200  and  m  =  400  grid  points. 
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Figure  12:  Ex.  6,  Woodward  and  Colella’s  double  blast  test.  The  density  profile  at  time  T  =  0.038  using  m  =  800 
(left)  and  m  =  1200  (right)  grid  points.  The  fine  resolution  reference  solution  is  shown  in  red. 


6.2.  Example  6:  Woodward  and  Colella’s  double  blast  test 

The  double  blast  test  of  Woodward  and  Colella  uses  the  initial  condition 

(pi  =  \,U]  =  0,pi  =  1000,  for*  <1/10, 
p,u,p  =  <pm  =  1  ,um  =  0 ,pm  =  1/100,  for  1/10  <  x  <  9/10, 

[pr  =  1,  ur  =  0 ,pr  =  100,  for  9/10  <  x. 

This  is  a  difficult  problem.  Nevertheless,  reasonably  good  results  are  obtained  by  our  CWEN03 
scheme  using  m  =  800  and  m  =  1200  grid  points,  as  shown  in  Fig.  12.  The  reference  solution  was 
computed  by  a  MUSCL  scheme  with  grid  size  m  =  4000. 

7.  Numerical  results  for  problems  in  two  space  dimensions 

We  now  discuss  some  numerical  examples  in  two  space  dimensions. 

7.1.  Example  7:  Rigid  body  rotation 

We  first  test  the  convergence  rate  of  the  scheme  as  applied  to  a  smooth  problem  involving  two 
dimensional  rigid  body  rotation  (cf.  [9]).  The  problem  is 

u,  -  ((y  -  1  )u)x  +  ((x  -  1  )u)y  =  0,  jc  e  [0, 2],  y  e  [0, 2],  t  >  0,  (57) 

with  a  smooth,  radial  bump  function  as  the  initial  condition,  defined  as 

u(x,y,  0)  =  |[«A(1  +  r(x,y))ifj{  \  -  r(x,y))  +  1], 
r{x,y)  =  y]{x  -  l)2  +  (y  -  l)2, 

where  ij/(s)  =  e~l,s2  for  5  >  0  and  if/(s)  =  0  otherwise. 

We  solve  the  problem  using  a  time  step  that  corresponds  to  a  fixed  CFL  number  of  0.1,  which 
is  At  =  0.1  h,  for  a  given  h  =  hx  =  hy  =  2/m  grid  elements.  In  Table  6  we  see  the  L*  and  L“  norms 
of  the  errors  and  the  corresponding  orders  of  convergence  as  m  increases.  Third  order  convergence 
is  observed,  as  expected. 
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Table  6:  Ex.  7,  smooth  rigid  body  rotation.  Errors  and  convergence  order  at  t  =  2n,  using  At  =  0. 1  h  (CFL=  0.1), 
where  h  =  hx  =  liy  =  2/m. 


m 

L'h  error 

order 

error 

order 

20 

1.62287E-02 

3.32619E-02 

40 

6.36774E-03 

1.34970 

1.58469E-02 

1.06967 

80 

1.04839E-03 

2.60260 

2.56955E-03 

2.62461 

160 

1.45170E-04 

2.85236 

1.83471E-04 

3.80790 

320 

1.80685E-05 

3.00620 

2.25254E-05 

3.02592 

640 

2.24239E-06 

3.01036 

2.72756E-06 

3.04587 

Figure  13:  Ex.  7,  rigid  body  rotation  of  a  square.  One  eighth  and  one  quarter  rotation  of  a  square  using  40  grid 
elements  in  each  coordinate  direction  and  At  is  0.425  of  the  CFL  limit. 


Figure  14:  Ex.  8,  a  Burgers’  problem  in  two-dimensions.  The  solution  at  time  T  =  1.5  is  shown  in  profile  and  as 
contours  using  It  -  1/80  for  the  two  plots  on  the  left  and  h  -  1  /320  on  the  right,  with  At  being  0.425  of  the  CFL  step. 


We  next  show  the  solution  for  a  rigid  body  rotation  of  a  square  (cf.  [7,  8]),  to  see  how  well 
the  scheme  handles  discontinuous  solutions.  The  domain  is  [0, 2]2,  and  the  solution  vanishes 
outside  the  smaller  square  [1/2, 3/2]2.  Using  40  grid  elements  in  each  coordinate  direction  and 
with  A t  being  0.425  of  the  CFL  limit,  the  solution  is  shown  in  Fig.  13  after  an  eight  and  a  quarter 
revolution.  The  discontinuity  is  captured  relatively  well,  with  no  spurious  oscillations,  even  though 
the  method  is  only  formally  third  order  accurate  and  the  mesh  is  relatively  coarse  (the  results  are 
comparable  to  those  in  [8]). 
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7.2.  Example  8:  A  Burgers’  problem  in  two-dimensions 

The  last  example  is  the  nonlinear,  two-dimensional  Burgers  equation 

du  d  / 1  7\  d  1 1  0\ 

~dt  +  ~dx\2U  j  +  dy\2U  )  =  x  6  ^  ?  >  (58) 

with  the  initial  condition  u(x,y,  0)  =  s i n 2 i /r.r)  sin2(7rv),  using  periodic  boundary  conditions.  In 
Fig.  14  we  present  the  solution  obtained  at  time  T  =  1.5  with  two  different  mesh  spacings,  80 
and  320  elements  in  each  coordinate  direction.  The  time  step  At  is  chosen  to  correspond  to  a  CFL 
number  of  0.425.  One  can  easily  see  that,  although  shocks  develop  in  the  solution,  they  are  well 
resolved  and  exhibit  no  spurious  oscillations.  The  result  is  comparable  to  [8,  Fig.  4.5],  although  in 
that  paper  it  appears  that  the  solution  has  been  plotted  in  reverse  in  both  directions  (and  they  show 
only  one-quarter  of  what  we  plot  here). 

8.  Conclusions 

We  have  presented  a  WENO  re-averaging  (or  re-mapping)  technique  that  converts  function 
averages  to  a  new  grid  and  to  high  order.  Nonlinear  weighting  ensures  that  we  maintain  the  es¬ 
sentially  non-oscillatory  property  of  the  re-averaged  function  values.  We  call  the  new  grid  the 
reconstruction  grid,  since  we  use  it  to  obtain  other,  standard  high  order  WENO  reconstructions 
of  the  function  averages,  such  as  high  order  point  values.  By  choosing  the  reconstruction  grid 
to  include  a  point  of  interest,  we  can  reconstruct  a  high  order  function  value  using  only  positive 
weights. 

We  applied  the  re-averaging  technique  to  define  a  classic  CWEN03  scheme  in  one  space  di¬ 
mension  that  combines  two  linear  polynomials  to  obtain  formal  third  order  accuracy.  Such  a 
scheme  cannot  otherwise  be  defined,  due  to  the  nonexistence  of  linear  weights  for  high  order  re¬ 
construction  at  the  center  of  a  grid  element.  Our  new  scheme  uses  a  compact  stencil  of  three 
solution  averages  in  the  reconstrutions  (and  five  for  the  overall  scheme,  including  the  Runge-Kutta 
step).  Stencils  are  combined  using  only  positive  linear  weights  in  both  the  re-averaging  step  and  in 
the  step  for  function  reconstruction  at  the  center  of  the  computational  grid  element  (which  is  now 
at  a  grid  point  of  the  reconstruction  grid).  There  are  actually  two  variants,  depending  on  whether 
one  chooses  to  use  a  uniform  or  a  nonuniform  reconstruction  grid. 

The  new  CWEN03  scheme  is  different  from  what  we  call  the  augmented  CWEN03  scheme 
of  Levy,  Puppo,  and  Russo  [7,  8],  since  their  scheme  uses  directly  the  quadratic  approximation  in 
addition  to  the  linear  ones. 

There  are  classic  CWENO  schemes  of  order  five,  nine,  etc.,  so  we  do  not  discuss  these  orders. 
Prior  to  this  work,  there  were  are  no  classic  CWENO  schemes  of  order  three,  seven,  eleven, 
etc.  due  to  the  nonexistence  of  the  linear  weights  for  reconstruction  at  the  grid  element  center 
points  for  these  orders  of  accuracy.  The  CWEN03  scheme  was  given  in  detail,  and  the  novel 
re-averaging  and  reconstruction  at  center  points  idea  that  we  presented  allows  one  to  define  the 
other  order  schemes  in  a  straightforward  way.  The  key  is  to  be  able  to  re-average  the  function  to 
a  reconstruction  grid  to  any  order  of  accuracy.  A  proof  that  this  can  indeed  be  done,  in  fact  on 
nonuniform  grids,  is  presented  in  [15,  10]. 
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We  extended  the  scheme  to  approximate  problems  with  two  space  dimensions.  The  extension 
was  developed  easily  in  an  essentially  tensor  product  form,  iterating  the  one  space  dimensional 
scheme  and  maintaining  the  compactness  of  the  stencil.  The  extension  is  therefore  relatively 
straightforward  to  implement,  and  it  is  computationally  efficient.  We  proved  that  the  scheme 
maintains  formal  third  order  accuracy.  Due  to  the  way  that  the  two-dimensional  scheme  is  defined, 
extension  to  problems  in  higher  space  dimensions  is  straightforward  to  implement  as  well. 

The  numerical  results  show  that  our  CWEN03  scheme  behaves  as  expected.  It  is  third  order 
accurate  for  smooth  problems.  It  also  gives  good  results  for  non-smooth  problems,  including 
Riemann  problems  for  the  Burgers,  Buckley-Leverett,  and  Euler  equations. 
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